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Abstract 

A new description of the binary fluid problem via the lattice Boltzmann method is 
presented which highlights the use of the moments in constructing two equilibrium 
distribution functions. This offers a number of benefits, including better isotropy, 
and a more natural route to the inclusion of multiple relaxation times for the binary 
fluid problem. In addition, the implementation of solid colloidal particles suspended 
in the binary mixture is addressed, which extends the solid-fluid boundary con- 
ditions for mass and momentum to include a single conserved compositional order 
parameter. A number of simple benchmark problems involving a single particle at or 
near a fluid-fluid interface are undertaken and show good agreement with available 
theoretical or numerical results. 
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1 Introduction 



The lattice Boltzmann equation (LBE) offers a very attractive way to study 
complex fluid flows governed by the Navier-Stokes equations (for a recent re- 
view see, e.g., [1]). This is particularly the case for flows involving irregular 
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and/or changing geometry, e.g., flows in porous media, multiphase systems, 
or multicomponent fluid mixtures. In the case of fluid mixtures where flows 
can have highly convoluted evolving structure, the LBE approach avoids the 
need to track the time evolution of the interfaces between two or more com- 
ponents — pinch-offs and topological reorganisations occur naturally on the 
lattice. The situation becomes more complex still if solid particles, or col- 
loids, are suspended in the flow. Such colloidal suspensions exhibit a range 
of interesting properties, and are of considerable interest in many technolog- 
ical and industrial applications (e.g., the food, cosmetic, and pharmaceutical 
industries) . 

A number of LBE approaches to multiphase and multicomponent fluids have 
been advanced [2,3], which adopt different strategies to capture the important 
thermodynamic interactions which give rise to surface tension between the 
phases. While numerical LBE results for single phase flow are very robust, 
multiphase flow is more problematic: if the thermodynamics is not consistent, 
spurious unphysical currents can be generated in interfacial regions owing to 
lack of detailed balance. While a complete, thermodynamically consistent, 
LBE approach to the binary fluid problem is still an open issue, the current 
methods can be used to provide useful results for many problems of interest, 
provided the errors incurred are monitored carefully, e.g., for phase separa- 
tion dynamics [4,5] and droplet motion and break-up in shear flow [6]. An 
alternative approach [7] is to use the LBE in the momentum sector, but em- 
ploy a standard finite difference technique for the evolution of a scalar order 
parameter describing the composition. 

Solid-fluid boundary conditions for stationary objects of complex geometry 
are readily implemented within the LBE. In addition, a general method for 
the inclusion of moving solid objects in a single phase is available [8] which 
allows the representation of, e.g., moving spheres [9] and eUipsoids [10]. These 
objects are generally resolved on the fluid lattice (they are at least several 
lattice spacings in size) without the need for complicated curved boundary 
treatments. Long range hydrodynamic interactions between particles are cap- 
tured by the LBE and, even when particles are separated by distances small 
compared to the lattice size, potentially important short-range lubrication in- 
teractions can be included for simple objects such as spheres by adding back 
any unresolved contribution using an analytical expression [11,12]. 

In this work, a number of these different strands are combined and extended to 
address a new problem: that of colloids in a binary solvent. First, a binary fluid 
LBE following [2] is adopted in which the equilibrium distributions are cast in 
a different, and more intuitive form. This leads to a number of improvements in 
the behaviour of fluid-only problems. Second, solid-fluid boundary conditions 
are extended to include colloidal particles. As one important advantage of the 
LB approach is the relative ease with which solid-fluid boundary conditions are 
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included, it is useful to retain LB in both the momentum and thermodynamic 
sectors. The boundary conditions for the two sectors can then be implemented 
in a consistent fashion. 



The paper is organised as follows. The following section gives an overview of 
the basic LBE algorithm used and the extension of solid-fluid boundary con- 
ditions for moving particles to the order parameter sector. Section 3 presents 
a number of results for simple test problems for the fluid only, and including 
a single particle. 



2 The Lattice Boltzmann Approach 



In this section, the LB approach to the binary fluid problem is described. As 

this differs from previous approaches [2,5], so the formalism is first set out 
for the lattice-Boltzmann equation (LBE) for the single phase in a way which 
makes clear the extension to the lattice kinetic equation (LKE) describing the 
second phase. 



2. 1 LBE for Navier-Stokes 

The LBE is commonly used to solve the conservation laws for mass and mo- 
mentum describing the hydrodynamics of an isothermal fluid 

dtp + V.g^O (1) 

and 

dtg + v.n = 0. (2) 

These equations describe the dynamics of the mass density p and the momen- 
tum density g = pv for fluid velocity v. The momentum density can be written 
as the mass flux ga and the momentum flux as Ua/s so that for a Newtonian 
fluid 

na/3 = QaVp + pSa/3 " vC^aVf) + V pV^) - C^aji^ pVp, (3) 

where Greek subscripts denote Cartesian components. Together, these equa- 
tions lead to the Navier-Stokes equations for isothermal flow, namely 

p(9tv + v.Vv) = -Vp + r]V\ + CV(V.v), (4) 
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where p is the pressure, while rj is the shear viscosity and ( is the bulk viscosity. 



The Lattice Boltzmann equation is derived from the (continuous time) discrete 
velocity equation 

dtfi + Ci.Vfi^-Cij{fj-fJ') (5) 



where the right-hand side corresponds to a linearised collision operator. The 
discrete velocity set {cj} are the nodes of a Gauss-Hermite quadrature [3,14] 
which ensure that the conserved moments of the distribution function, /j, have 
the same values as in the continuum. For the isothermal LBE, these are the 
mass and momentum densities 

P = E/^> (6) 

i 

Qa^^fiCia- (7) 

i 

The Boltzmann kinetic description is restricted to a dilute gas with an ideal 
equation of state p = n/c^T = pc^, where n is the number density and 
is the isothermal speed of sound. It is convenient to subtract the trivial ki- 
netic contribution to the pressure from the momentum flux tensor to define a 
deviatoric momentum flux 



where Uaf) = Ej fiCiaCip is the momentum flux and Qia/s = CiaCip - c^Sap is 
referred to as the kinetic projector [1]. 

In the most commonly used DdQn models, with n velocities (or quadrature 
nodes) in d dimensions, the equilibrium distribution functions are given by 



fi -Wi 



, pVaCia , PVaV/sQiafS 
P+ ^ + 



2cf 



(9) 



where the Wi are weights defining the quadrature and repeated Greek indices 
are understood to be summed over. This form is obtained by a truncation 
of the Hermite polynomial expansion of the Maxwell- Boltzmann distribution 
at second order [14]. Finally, the BGK coUision matrix Cij must satisfy the 
constraints imposed by the conservation laws and rotational symmetry. In the 
single relaxation time approximation Cij = Sij/r, and the shear and bulk 
viscosities are related to the relaxation time t hj rj = pc^r and ( = (2/(i)pc^T, 
respectively. For multiple relaxation time models, the above relations remain 
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valid if T is replaced by separate relaxation time for the shear and bulk viscous 
stress. 

The final LBGK equation is obtained by a further discretisation of the contin- 
uous time equation. A second-order characteristic-based method [15] can be 
followed to obtain 

/,(r + cAt; t + At)- Mr- 1) = -(/,- /r)/(r + f ). (10) 

This does not introduce lattice artefacts at second order so that the definition 
of the viscosity remains as in the continuum, unlike that in a first order scheme, 
where the lattice error is absorbed into the viscosity to give rj — pc1{t — 1/2). 



2.2 LKE for Cahn-Hilliard 



The starting point for the binary fluid mixture in this work is the (Landau- 
Ginzburg- Wilson) free energy functional which describes the total energy of 
a system of fixed volume as a functional of a single compositional order 
parameter (f){r,t). The order parameter measures the ratio of the number 
density of particles of the two components of the mixture rii and n2 so that 
= (ni — n2)/(ni -|- ^2). The choice of the free energy uniquely defines the 
physical quantities of interest such as the fiuid-fiuid interfacial tension a and 
the interfacial width ^o- The chemical potential is obtained from the free energy 
functional via 



so that the thermodynamic force density acting on the fiuid is then — 0V//. 
The equation of motion for the order parameter is the Cahn-Hilliard equation 
dt(l) + V.((j)V - MVii) = 0. (12) 



This is a conservation law involving the single conserved quantity 0, the flux 
of which 0v — MVyU is made up of an advective component related to the 
fluid velocity v and a diffusive component related to the gradient of the chem- 
ical potential by an order parameter mobility M. The mobility is assumed 
to be constant and independent of 0. A kinetic relaxation scheme for the 
Cahn-Hilliard equation is obtained by introducing distributions Qi obeying 
the kinetic relaxation equation 

dtQi + Ci.Vgi = -Cfj {gj - gf). (13) 
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The discrete velocities Cj are the same as for the LBE, and physical quan- 
tities are again related to moments of the distribution via = ^iQi and 
jo, — J2i9iCia- In choosing the form of the equilibrium distribution, constraints 
similar to those of [2] are employed, namely 



i 

i 



(14) 
(15) 
(16) 



While the physical interpretation of the constraints on the zeroth and first 
moments is clear, that for the second moment is somewhat less so. There are 
two contributions, the first of which ensures a diffusive contribution to the 
evolution of the order parameter related to the chemical potential, while the 
second represents an advective flux related to the velocity held. It is then 
possible to write down an equilibrium distribution by analogy with Eq. 9, i.e.. 



-.eq 



Wi 



0H 5- H 



(17) 



While this choice satisfies the above constraints, it is not unique. As it turns 
out, a slightly modified version is required in practice (see following section). 

Further discretisation again provides an LBGK equation with a single relax- 
ation time for the order parameter 

gi{v + c,At; t + At) - g,{r- 1) = -{g, - g^Mir'^' + f ). (18) 



The order parameter mobility is related to the relaxation time via M = c'^t'^, 
which again can be adjusted to absorb the discrete lattice correction. Note 
that the relaxation time is not fixed in this approach, in contrast to previous 
work. 



2.3 Implementation 

The choice of free energy functional follows that of [5] , where 

F[(/>] = / dr[| + lBcf>' + !«:( V0)2] . (19) 
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This describes a symmetric binary mixture in which the bulk free energy is 
related to two parameters A and B, while a term in a third parameter k 
penalises gradients in the order parameter, i.e., it tends to minimise curvature 
of the interface. The fluid-fluid interfacial tension is then 

a = (-8«;AV9S2)V2. (20) 



If the further constraint that that ^4 = — S < is added, then the order pa- 
rameter lies predominantly on the interval [—1,1]. The equilibrium interfacial 
profile is a tank with characteristic width 

C = (2«;/A)V2. (21) 



For the results presented in this work, the thermodynamic force arising from 
the order parameter sector is introduced via a correction to the equilibrium 
stress from a chemical pressure tensor 



a/3 



(22) 



Having made these choices, the lattice Boltzmann algorithm of collision fol- 
lowed by propagation can be implemented in the normal way. The post- 
collision distributions, denoted /*, are based on Eq. 9 so that 



ft 



+ 



(23) 



where the density and momentum are unchanged by the collision process so 
that p* = p and ga = Qa- The deviatoric stress is relaxed with a single relax- 
ation time satisfying 77 = pc^r 

Sip = Sap - {Sap - Sll)/{T + f ) (24) 



where = pUaUp+Pap- Note that this reprojection of the physical properties 
to the distribution readily admits the use of multiple relaxation times, e.g., if 
separate values for the shear and bulk viscosities are required. 

For the order parameter, there is only one conserved quantity in the collision 
(f)* — (f), where again the star refers to post-coUision quantities. The order 
parameter flux is relaxed with single relaxation time satisfying M — c^t'^ so 
that 

j:-ja-ija-fJ)/{rUf) (25) 
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with j^^ — (f)Va- Finally, the values of and v arc known, allowing a repro- 
jection of the physical quantities to the g* distributions. If the equilibrium 
distributions Eq. 17 are used the model is, in general, numerically unstable. A 
simple finite-difference expansion [13] of the evolution equation (not shown) 
shows a spurious up-gradient diffusion arises which is the cause of the insta- 
bility. Instead, it is possible to use 



where the 6io has the effect of moving most of the order parameter into the 
non-propagating rest distribution qq. This choice is found to have very good 
stability properties, while also satisfying the constraints on the moments of Qi. 



2.4 Colloidal particles 

A very general method for the representation of solid particles within the LB 
approach was first put forward by Ladd [8]. Solid particles (of any shape) can 
be defined by a boundary surface which intersects a set of vectors {cijAt} which 
connect lattice nodes inside and outside the surface. These are referred to as 
boundary links. In the original approach, this set of links defined a (spherical) 
shell with fluid occupying all lattice nodes both inside and outside the "solid" 
object. While this so-called internal fluid can be criticised as unphysical, it 
actually exerts significant effect on the dynamics of the particle in a single 
phase flow only on the time scale that it takes sound waves to cross the particle 
[17]. However, in a binary fluid mixture, it becomes essential to have truly 
solid particles to prevent possible unphysical transfer of fluid across particle 
surfaces. Such unphysical transfer would impact on the net composition of 
the real fluid (i.e., that outside the particles). This is particularly true for 
particles which wet one species of fluid preferentially (non-neutral wetting), 
where a particle might capture one species of fluid and hence unbalance the 
net composition and/or generate spurious thermodynamic forces in the region 
of the surface. A number of implementations of truly solid particles have been 
developed [16,17]; this work extends that of Nguyen and Ladd [11] to the 
binary fluid problem. 

A schematic diagram showing the distribution of links for a section of a spher- 
ical particle is shown in Figure 1. The centre of a sphere of radius a at Fq 

deflnes the position of the links, and the particle moves smoothly across the 
lattice with linear velocity U and angular velocity ft. Boundary nodes are 
deflned to be half way along the links, which the set of vectors joining the 
centre of the sphere to the boundary nodes denoted {vb}. For a single integral 



9* = 4>*Sio + Wi 




+ 



2cl 



(26) 
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Fig. 1. A schematic showing the distribution of boundary hnks for a section of a 
spherical particle. The nominal radius of the particle is a, while the vector joining 
the centre to the boundary nodes is r^. If the particle moves relative to the lattice, 
fluid nodes may be exposed (right), in which case fluid with appropriate properties 
must be added back. The actual lattice used in calculations is D3Q15. 

particle (having a full complement of links) 

J2^c,c, = (27) 

b 



and 



b 



w. 



jb X Cfc) = 



(28) 



where Wc^ are again the quadrature weights appropriate for the boundary links. 
These results will be useful in the next section. 



2.4-1 Boundary conditions and particle dynamics 

The boundary condition developed by Ladd is usually referred to as bounce- 
back on links (BBL). For a moving particle, this takes the incoming distribu- 
tion at the solid-fluid boundary, ff,, and reflects it back along the incoming 
direction 

fb' = fb- 2wc^pub.Cb/c1 (29) 



where Cb' = — Cf,. The correction is related to the local solid body velocity at 
the boundary node 

Ub = V + nxvb, (30) 



and has the effect of transferring mass from the leading edge of the particle 
to the trailing edge. If the approximation is made that p ^ pq, i.e., that the 
density is approximately constant around the particle, then Equations 27-30 
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can be combined to show that the BBL apphed around the whole particle 
conserves mass. 

Conservation of linear and angular momentum entails computing the momen- 
tum transfer from the fluid to the particle as a result of BBL, and updating the 

particle Uncar and angular velocity appropriately. For this, the fully implicit 
method described by Nguyen and Ladd [11] is used. 

The BBL condition for density is extended to the order parameter [18] so that, 
at the surface, the distributions gb are again reflected with a correction 

Qb' = gb- 2Wc,0U(,.Cb/c^. (31) 

However, the order parameter around the particle cannot, in general, be 
approximated by a constant as is the case for density. This means that the 
order parameter is not conserved at BBL by an amount 

50 = 2 ^ Wc,(j)Ub.Cb/cl (32) 

b 

for a given particle. To ensure the composition of the fluid does not drift over 
time, the deficit or excess 50 is added back as a correction at the following 
time step: 

Qb' ^ Qb- 2wc(,0U6.Cb/c^ - Wc^5(j)/ '^cb- (33) 

b 

The size of the additional term is generally small compared with the term in 
Ufe.Cfe, so this represents a negligible perturbation to the motion. 

Perhaps more serious are the corrections that arise when the particle changes 
shape as it moves across the lattice. When the particle motion exposes a fluid 
lattice node or recovers a solid node, fluid with the appropriate properties 
must be added or removed. For density, this gives rise to a correction to the 
BBL Eq. 29 that ensures the mean fluid density po does not drift [11]. A 
similar correction can be made in the order parameter sector which ensures 
that the mean fluid composition 0o does not change. For example, if fluid is 
added at a newly exposed lattice node, then new properties are determined 
by a linear interpolation of the distribution functions at adjacent fluid sites. 
For order parameter, an extra correction to the BBL at the following step of 
S(f) = (0 — 0o) is required to maintain mean fluid composition, being the 
order parameter added or removed. These corrections lead to small unphysical 
fluctuations in the order parameter near a moving colloid surface. At present, 
these are accepted as the cost of preventing drift in the fluid composition. 
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Finally, the gradient in the order parameter field is required to compute Pa/3 
from Eq. 22. Near the solid particles, these gradients are computed following 
[18] where the order parameter is extrapolated along links to the boundary 
nodes. In this work only neutral wetting is considered where the contact angle 
between the sohd and the fluid-fluid interface is 90 degrees. It should be noted 
that in imparting the thermodynamic force to the fluid via the stress Pap, 
there is no direct thermodynamic force on the colloid. The order parameter 
only then affects the particle motion indirectly via the fluid velocity fluid. 

2.4-2 Particles close to contact 

The BBL for the density distribution functions fi allows the net hydrody- 
namic force and torque on a particle to be computed, from which the particle 
velocity and position can be updated in turn. For an isolated particle this is 
straightforward. However, if two or more moving particles are close enough on 
the lattice scale that there are no fluid sites in the interstice, the full lubrica- 
tion force between the particles will not be resolved. In this case, it is possible 
to add back the unresolved part from the analytical expressions available for 
pairwise forces between spheres [21] and, after appropriate calibration, recover 
the correct lubricating behaviour at close approach [11]. 

The same effect can occur in the order parameter sector when, for example, 
two particles are close together at a fluid-fluid interface. The particles should 
experience capillary forces owing to the curvature of the interface in the gap 
between them, but, again, if no fluid is present this force will not be resolved. 
Unfortunately, there is no simple way to add back the unresolved force as 
is done for the hydro dynamic lubrication as the capillary force depends in a 
complex way on the curvature of the interface in the vicinity of the particles. 
While these capillary forces may be important for some problems, they do not 
have the same impact as the lubrication forces in, for example, preventing the 
particles overlapping. However, the results presenting in the following section 
only consider a single particle. Calibration of the hydrodynamic radius of 
different particles for the single fluid is carried out following [11]; these are 
assumed to remain unchanged for the binary fluid. 



3 Results 

3.1 Fluid only 

As a demonstration that the equilibrium distributions presented in the pre- 
vious section lead to better behaviour in the fluid sector than those used 
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Fig. 2. An initially steady spherical interface between two fluids is allowed to relax 
for some 100s of time steps. The LBE method described in [2] ("SOY") is compared 
with the new method. While the two methods give similar results at surface tension 
0.0042(a and b), anisotropics are exposed using the SOY distributions at higher 
surface tension 0.055 (c and d). 

previously, Fig. 2 shows a comparison of results for a simple test problem. An 
initially steady spherical fluid droplet of one phase is initialised surrounded by 
the second phase. The system, here 64^ lattice sites with periodic boundary 
conditions, is then allowed to relax for a few hundred time steps. This problem 
is repeated for two parameter sets taken from Table 3 of [5] corresponding to 
surface tensions of cr = 0.0042 and a = 0.055. While there is little visible 
difference between the results for the lower surface tension, anisotropy in the 
old distributions is exposed at higher surface tensions which do not infect the 
new distributions. This improved behaviour is also manifest in improved nu- 
merical stability. While by no means unconditionally stable, use of the new 
distributions allows a wider range of parameter space to be investigated safely. 



3.2 Single particle at fluid-fluid interface 



A spherical particle subject to gravity may be suspended at a fluid-fluid in- 
terface by interfacial tension. For a particle which has a density difference Ap 
with the surrounding fluid, the key dimensionless number in this situation is 
the Bond number 

Bo = a'^Apg/a (34) 



where g is the gravitational acceleration. The Bond number expresses the bal- 
ance between the downward (or upward for buoyant particles) force and the 
opposing interfacial tension. If Bo << 0(1) the interfacial tension can sup- 
port the particle with little deformation. As Bo increases, the interface bows 
downward until the particle can no longer be supported. For a neutrally wet- 
ting particle with a contact angle of 90° between solid and fluid-fluid interface, 
the critical Bond number is 3/4; particles break away from the interface for 
higher values. 

In the LB model, a single particle of radius a is placed at rest across an initially 
flat interface and subject to a downward force. The particle is allowed to fall 



12 



n3 

N 



0.0 



-0.5 



-1.0 



-1.5 
10 



10 10 10 

Bond Number 



10° 



Fig. 3. The normalised equilibrium displacement zja of a single particle below an 
initially flat interface as a function of Bond number. The different symbols represent 
particles of different radius: a = 2.3 circles, a = 3.71 triangles, and a = 4.77 squares. 
The curve is the theoretical result [19] is valid for Bo << 1. The particle detaches 
from the interface for Bond numbers higher than the critical value B = 3/4. 

until it comes to rest, in which case the equilibrium displacement below the 
level interface is measured, or it detaches from the interface. The width of the 
system in the horizontal direction is at least 10a so that there is no significant 
curvature of the interface at the periodic boundary conditions. In addition, the 
second interface required in the system is placed at sufficient distance that it 
does not affect the motion of the particle. 

Figure 3 shows the normalised displacement /i/a as a function of Bond number 
for a number of different particle sizes. For Bo < 0.01 the interface is essen- 
tially rigid with no significant displacement, while the displacement increases 
to a significant fraction of the particle size for higher Bond numbers. The 
agreement between the numerical results and the theoretical result is gener- 
ally excellent. As Bo approaches the critical value of 3/4 the smaller particles 
drop somewhat below the theoretical curve and can detach from the inter- 
face at Bond numbers as low as 0.55. This is likely to be caused by a poor 
representation of the contact line for smaller discrete particles. 



3.3 Particle approaching a fluid- fluid interface 



As a second test of particle motion in a binary mixture, the drag on a sin- 
gle sphere sedimenting toward a stationary fluid-fluid interface is computed 
(Fig. 4). First, the mean drag on the sphere is obtained in a single phase LB 
calculation to provide the calibration Girrja as for the hydrodynamic radius 
(the nominal radius used here is a = 2.3). The drag is then recomputed for 
the same particle in the binary fluid as it sediments vertically toward a hori- 
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Fig. 4. The drag (normalised by Gvrr/a) on a single sphere sedimenting toward a 
stationary fluid- fluid interface as a function of the normalised separation h/a. The 
solid line is based on tabulated data from [20], while the symbols represent the 
current model: triangles represent interfacial width ^^"^ = 0.8 and crosses ^ = 1.6 
lattice units. 

zontal fluid-fluid interface. The Bond number is kept small enough that there 
is little deformation in the interface as the particle approaches. The drag on 
the sphere is measured as a function of the normal separation h/a; as the ex- 
act value depends upon the position of the particle relative to the lattice, the 
symbols in Fig. 4 represent an average over 20 different particle configurations 
and binned at intervals in h/a. 

The results show that, far from the interface, the drag is very similar to that 
seen in the single phase. This suggests the redistribution of order parameter 
as the particle moves does not strongly influence the dynamics. As the parti- 
cle approaches the interface, the drag increases but then decreases sharply as 
the leading edge of the particle touches the interface. The effect of differing 
interfacial width is shown: the sharper interface induces a steeper rise in the 
drag as the particle approaches. However, in both cases the particle is flnally 
captured by the interface. For comparison, the exact results for a similar prob- 
lem in which a sphere approaches a flat interface of zero thickness is shown 
[20]. In this case, the particle never actually reaches the interface. 



4 Closing Remarks 

This work has demonstrated the use of a new binary fluid lattice Boltzmann 
approach applied to the problem of colloidal particles The representation of 
colloidal particles within an LBE binary fluid allows a large number of inter- 
esting physical problems to be investigated. The results presented shown an 
excellent agreement with exact results for a number of simple test problems, 
and provides confldence that the approximations made in the approach do not 
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result in undue errors. 



A number of potential improvements can be identified. First, the thermody- 
namic force follows [2] in entering via the equilibrium stress. This can be 
added directly as an additional force on the fluid in the LBE. The implemen- 
tation of this additional force in the presence of solid particles requires some 
care and is addressed elsewhere [22]. Second, the problem of resolving near- 
contact capillary forces would be most elegantly addressed via the use of a 
finite-volume like approach where there is always some fiuid retained in the 
interstice between particles. This approach would also have the benefit of pre- 
venting abrupt changes in the discrete shape of a moving particle. However, 
this would represent a considerable increase in complexity over the current 
approach. 

In the meantime, the work demonstrates the flexibility of the lattice Boltz- 
mann method in addressing problems that are very difficult to address at all, 
and currently intractable using other methods. 
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